Foam density mapping via THz imaging

Plastic foams, near-ubiquitous in everyday life and industry, show properties that depend primarily on density. Density measurement, although straightforward in principle, is not always easy. As such, while several methods are available, plastic foam industry is not yet supported with a standard technique that effectively enables to control density maps. To overcome this issue, this paper proposes Terahertz (THz) time-of-flight imaging using normal reflection measurements as a fast, relatively cheap, contactless, non-destructive and non-dangerous way to map plastic foam density, based on the expected relationship between density and refractive index. The approach is demonstrated in the case of polypropylene foams. First, the relationship between the estimated effective refractive index and the polypropylene foam density is derived by characterizing a set of carefully crafted samples having uniform density in the range 70–900 kg/m3. The obtained calibration curve subtends a linear relationship between the density and the refractive index in the range of interest. This relationship is validated against a set of test samples, whose estimated average densities are consistent with the nominal ones, with an absolute error lower than 10 kg/m3 and a percentage error on the estimate of 5%. Exploiting the calibration curve, it is possible to build quantitative images depicting the spatial distribution of the sample density. THz images are able to reveal the non-uniform density distribution of some samples, which cannot be appreciated from visual inspection. Finally, the complex spatial density pattern of a graded foam sample is characterized and quantitatively compared with the density map obtained via X-ray microscopy. The comparison confirms that the proposed THz approach successfully determines the density pattern with an accuracy and a spatial scale variability compliant with those commonly required for plastic foam density estimate.

compared to µ-CT, N-CT is based on the enhanced neutron attenuation of polymers (due to the high hydrogen content in their macromolecules) and therefore the contrast for such structures is in general higher 10 .However, large instruments and complex post-processing efforts are needed also in this case 11 .To our knowledge, it seems that the plastic foam industry is not supported with an effective technique that enables a routinely determination and/or control of the density maps.
THz are electromagnetic waves ranging from 0.1 to 10 THz (wavelength from 3 mm to 30 μm), which can penetrate a wide range of non-conducting materials, e.g., plastics, polymers, ceramics, wood, and glass 12 .Accordingly, they allow for a contactless, safe, and non-destructive characterization of these materials, not limited to the surface of the object under test but also to its interior.Typically, structures of thickness up to a few centimeters along the incidence direction can be investigated, providing information at millimeter scale.THz imaging and spectroscopy have been considered in many fields [13][14][15][16][17] , among which cultural heritage 18 , security 19 and pharmaceutical applications 20,21 , only to cite a few examples.Plastic foams have been, also, characterized by using THz and the performed spectroscopy studies have shown that their properties change depending on the kind of foam, blowing agent and density 22 .As a novel contribute to the available literature, this paper investigates the use of THz imaging to obtain quantitative density maps of graded plastic foams, which are attracting a growing interest in the foam industry thanks to their peculiar mechanical properties directly depending on the designed 3D density map.
In the ever-growing body of literature, some studies considering the use of THz for density estimation are worth being mentioned.An early attempt to map density was pursued by Koch et al. 23 , who exploited transmission measurements and derived a calibration curve to relate the density of wood, as obtained by gravimetricvolumetric methods, with the THz absorption as measured for the known sample thickness.Wood density characterization via THz transmission measurements was also considered by Kashima et al. 24 , where a linear regression approach was proposed to predict density and moisture content from refractive indices and absorption coefficients obtained from THz transmission measurements.THz has been also considered for density and porosity measurements of pharmaceutical tablets 21,25 , as well as for imaging foam density and defects 26 .This latter paper uses reflection measurements and an external reference structure.However, to the best of our knowledge, THz imaging of foam specimens with non-uniform volumetric density distribution has never been explored before.
The approach we propose is based on time-of-flight (ToF) THz imaging.With respect to the literature 16,19,[21][22][23][24][25][26][27][28] , which mainly consider (single) transmission measurements, the adopted approach requires that the sample under test is positioned on a flat metallic plane and scanned along the x-y plane in normal reflection mode (sometimes referred to as double transmission).For each measurement point, the reflected signal is processed to determine the effective refractive index of the sample along the signal path.Specifically, the estimated effective refractive index is an average parameter depending on the sample dielectric features along the signal propagation path (i.e., the z-axis in our configuration) and its dispersive behaviour in the frequency range of the THz probing pulse.As such, the estimated equivalent refractive index accounts for the density variations occurring across the sample below the surface and along its depth.Then, provided the relationship between effective refractive index and mass density is available, the approach delivers a 2D map that represents the spatial distribution of the density, in which each pixel encodes the average density along the third dimension.Although it is very simple, the proposed approach is an effective and widely applicable procedure to estimate the refractive index, and thus the density, of uniform or not specimens, without requiring a priori information about their point by point thickness, while avoiding optimization procedures.Effective processing tools are available to estimate refractive index and thickness simultaneously 27,28 but they account for uniform sample and a (single) transmission measurement setup.
The approach is initially applied to characterize a set of uniform samples with variable nominal density values in order to establish the relationship between the effective refractive index and the density.The experimentally achieved calibration curve reveals a linear relationship between these two quantities within the range of interest.This relationship is in good agreement with that derived by using the medium theory given by Scheller et al. 29 and it is also validated by using another set of samples, referred to as test samples.THz surveys of the test sample reveal that their measured average density agrees with the one predicted by the calibration curve.In addition, THz images show that some of these samples have a not-uniform refractive index distribution.Hence, by means of the calibration curve, a quantitative description of their non-uniform density pattern, invisible to visual inspection, is obtained.Finally, the approach is exploited to image the density pattern of a foam sample having a 3D optimized density map for the three-point bending load condition 6 .The density map achieved via THz is validated with the one resulting from X-Ray microscopy.To the best of our knowledge, a similar comparison has never been presented before and the obtained qualitative and quantitative coherence supports the usefulness of the proposed approach as a simple strategy to exploit THz technology for mapping volumetric density variations, thus answering an open issue in foam industry.

THz time-of-flight imaging
THz time-of-flight (ToF) reflection imaging allows the collection of a large amount of information about the sample under test.At each measurement point, the object is probed by a pulsed signal and the reflected waveform is collected within a certain observation time window.The reflected waveform exhibits only one peak due to the air-sample interface in the case of a not penetrable sample; otherwise, it accounts for the sample surface and its inner electromagnetic discontinuities.Specifically, if the sample is homogeneous and sufficiently thin to be entirely penetrated, the reflected waveform exhibits two peaks corresponding to the beginning (upper side) and the end (bottom side) of the sample along the wave propagation path, see Fig. 1.The time delay between these two peaks depends on the thickness of the sample and its electromagnetic properties (i.e., its refractive index).
If the sample has a layered structure made up of electromagnetic different materials, the collected waveform appears as a pulse train where the temporal delay between successive pulses provides information about the sample stratigraphy.Therefore, a 3D image representing surface and interior of the sample can be obtained by scanning the whole sample 30 .
The ToF t is related to the distance d between the THz probes and the detected discontinuities as: v = c/n being the electromagnetic wave propagation velocity into the medium wherein the propagation occurs, c the light propagation velocity and n the refractive index of the medium along the wave path.
In our set-up, the sample is positioned on a flat metallic plane and probed in normal reflection mode, i.e., the incidence angle is zero, see Fig. 1.Let t M and t SM denote the ToF referred to the metallic plane reflection when the signal propagates in air and when it crosses the sample, respectively, and let t S be the ToF related to the reflection due to the air-sample interface.
For each point (x p , y p ) belonging to the sample, the thickness, �(x p , y p ) is readily estimated from Eq. (1), as: where �t SM = t M − t S (x p , y p ) is the time-delay associated to the propagation that occurs in air (i.e., n = 1 ).Hence, the map of depicts the sample thickness, appraising variations whose size is equal or higher than the spatial range resolution of the adopted THz system.Assuming an effective medium approximation for the sample under test, according to which the medium is approximated as a homogeneous and non-dispersive medium along the z direction.Each measurement point (x p , y p ) scanning the sample surface can be characterized by an effective refractive index, say n eff (x p , y p ) , obtained from Eqs. (1) and (2) as: Accordingly, the plot of n eff provides a 2D quantitative map of the refractive index along the scanning surface with a spatial resolution depending on the measurement offset, which takes into account point-by-point variations of the sample thickness.

Samples preparation
Among possible plastic foams, this study has been focused on Polypropylene (PP) foams.Two types of PP foams samples were produced: • a set of cylindrical samples having a nominally uniform density in the range 70-900 kg/m 3 ; • a graded-density sample, a beam with a 3D density map optimally shaped for three point bending.
The foam production process is schematized in Fig. 2a.
The cylindrical samples were obtained from neat PP preforms (disks), punch cut from sheets.The sheets were produced from the granules using a compression molding machine (model P300P, Collin Gmbh, Ebersberg, Germany) as follows.The mould was preheated at 200 °C and the granules as received were placed on a flat press table for 2 min without pressure on a polytetrafluoroethylene sheet.The pressure was then increased to 5 bar in 1 min, then to 10 bar in 1 min and to 20 bar in 1 min.After cooling (to 25 °C), the pressure is brought to ambient condition and the produced PP sheet was extracted.
The graded-density sample was produced as detailed in Iaccarino et al. 6 .Briefly, preforms were produced by assembling 1 mm-thick slabs cut from PP sheet, compression molded as previously described.Polypropylene (PP) was kindly supplied by Sinopec Zhenhai Refning & Chemical Company (Zhejiang Province, China) as granules (commercial product identified as E02ES).CO 2 was used as blowing agent and was provided by SOL S.p.A. (Monza, Italy).
(1) For the production of foamed cylinders, disks were placed in a three-piece mould with cylindrical cavities of 8.3 mm in height and 23.2 mm in diameter (Fig. 2c).The moulds were placed in the pressure vessel and a temperature of 160 °C was reached under vacuum to melt the polymer.The temperature was then increased to 146.5 °C, where saturation with CO 2 was carried out at 120 bar in 10 min (sorption time).Finally, the ball valve was opened, leading to foaming of thobtainedks in the mould, where the foam could be shaped to the cylindrical cavity and set by cooling.Using preforms of various initial sizes and volumes, foamed samples of various densities were obtained for further testing (Fig. 2d).
For the production of graded-density beam, the optimized preform was placed in a box-shaped mould (9 × 140 × 15 mm 3 ) and subjected to the same foaming procedure described above, but the sorption time was 30 s, instead of 10 min.

Density measurements and volumetric fractions calculation
The average density of the produced foams ( ρ f ) was obtained by using the following equation: where m a is the weight of the foam in air, m w is the weight of the foam in water, ρ a is the density of the air and ρ w is the density of the water.At room temperature ( 23 ± 2 • C), ρ a = 0.00120 ± 0.00002 g/cm 3 and ρ w = 0.9975 ± 0.0005 g/cm 332 .m a ( ± 0.0005 g) and m w ( ± 0.0005 g) were measured using an analytical balance (MS Semi-Micro Model, Mettler Toledo, Milan, Italy).The same procedure was adopted to measure the density of neat PP (i.e., unfoamed PP, C6), denoted as ρ p .
Then, the volumetric fraction of the polymer, χ p , can be computed as follows: The absolute error on the measured foam densities, δ ρ f , and the absolute error on the calculated volumetric fractions of the polymer, δ χ p , were obtained as described in the "Appendix".The parameters characterizing the cylindrical shaped samples, referred as C1-C7 and T1-T6 are reported in Table1.

Measurement system
THz data have been gathered by means of a custom time-domain system developed by MenloSystems 33 .This device generates and gathers broadband THz pulses and enables measurement in both transmission and reflection modes.Figure 3 shows the system arranged for the reflection measurements mode exploited in this study.
The core of the system is the optical sampling engine (OSE) module, which exploits asynchronous optical sampling technique to perform high-speed scanning over some nanoseconds of time delay without using a mechanical delay line 34 , pushing the spectral resolution to the region of hundreds of MHz 35 .The OSE module includes two pulsed femtosecond lasers where one ultrafast laser delivers the pump pulse for the emitter antenna and the other laser delivers the probe pulse for the detector antenna.The lasers operate at a locked repetition rate with a tunable difference, providing the optical pulses for THz emission and measurement.The laser pulses are delivered via optical fiber to the THz photoconductive antennas (PCA), specifically the TERA15-TX-FC and TERA15-RX-FC fiber-coupled photoconductive antennas are used.The PCA are inserted in a compact THz reflection head with integrated optics for high-performance measurements, which is mounted on imaging stages equipped with stepper motors.These motors enable movement along the x-axis and the y-axis, resulting in a scan area of 300 mm × 300 mm and a minimum measurement spatial offset of 120 µ m.The object under test has been positioned at the focal distance, i.e. at about 3 mm from the reflection head for the case at hand.
At each measurement point, the system collects signals within a 100 ps observation time window, which can be moved along a time scan range of 1 ns.The maximum depth that can be investigated for a 100 ps observation time window is d max < 100 • 10 −12 • v/2 .This means that in the most favorable cases (c ∼ v) d max is less than 1.5 cm.Although the nominal bandwidth of the system is about 4 THz, in the case of normal reflection measurement and uncontrolled environmental conditions, the actual frequency range is B = [0.2− 2] THz.This means that, when the propagation occurs in air, the expected range resolution is r = c 2B = 83.3µm.

Measurement procedure
The measurements procedure is summarized in the following steps: 1.The samples are placed on a metal plate at the focal distance from THz reflection head.The scanning area is 30 mm ×30 mm for samples C1-C7 and T1-T6 and 100 mm ×30 mm for the graded-density sample.The spatial offset along x and y axis is 0.5 mm for C1-C7 and T1-T6 and 0.25 mm for the graded-density sample.At each measurement point, the stored waveform is the average of those measured in a time window of 1 s. 2. The raw data are processed by means of an ad-hoc filtering strategy 36 .First, a frequency band-pass filter (BPF) procedure is applied to remove low-and high-frequency signal components, by picking signal components in the actual [0.2-2] THz range.Then, a noise filtering procedure, based on the singular value decomposi-  tion (SVD) of the band pass filtered data matrix, filters out further noise components, without affecting the spatial resolution.3. The ToF t M is estimated by considering a generic processed waveform referred to a measurement point external to the sample and by looking for the time position of the maximum amplitude of the waveform, i.e. the peak due to the reflection by the metal (note that the metal is flat and parallel to the measurement aperture); 4. The ToFs t S and t SM are estimated for each measurement point (x p , y p ) belonging to the sample by searching the time positions at which the processed waveform exhibits its maximum amplitude values before and after t M .The maximum amplitude occurring before t M represents the peak due to the air-sample interface and its time position is t S , the other maximum amplitude is the peak due to the sample-metal interface and its time position is t SM . 5. The ToFs t S , t SM , and t M are used to compute the thickness value �(x p , y p ) by means of Eq. ( 2) and the effec- tive refractive index n eff through Eq. ( 3).

Calibration curve and THz density maps
To turn the effective refractive index into the density map the relationship between these quantities is experimentally established by analyzing the uniform density samples C1-C7, which show different and known density values, ranging from very low values up to that of neat PP.
To verify the measurement repeatability and the possible dependence of the results on the probed sample side, the samples have been analyzed by surveying both the largest sides and two data sets were collected for each side.Hence, for each sample, 4 effective refractive index and thickness maps were obtained, � i (x p , y p ) and n eff ,i (x p , y p ) with i = 1, .., 4.
The average effective refractive index n avg was estimated as follows.For the i-th dataset, n eff ,i (x p , y p ) was averaged among all measurement points to obtain the effective refractive index characterizing the dataset.Then, n avg was obtained as the average among the four values and associated to the density value ρ f corresponding to the sample at hand, as taken from Table 1.The average thickness avg was computed using the same procedure.
The calibration curve was then obtained performing a linear fitting of the average effective refractive index n avg estimated for all samples C1-C7, taking into account the standard deviation σ n arising from the computa- tion of n avg for each sample.
The calibration curve is coherent with the relationship predicted by the medium theory of Scheller et al. 29 .Moreover, it has been validated against the T1-T6 samples, to check the matching between the density as estimated from their average effective refractive index and the one measured in "Density measurements and volumetric fractions calculation" section.Also for these samples, four data sets are collected (two for each side).
For the T1-T6 samples, the density maps obtained from the effective refractive index by using the calibration curve are inspected to show the capability of THz ToF imaging of characterizing the spatial variability of the foam density.

Validation on the graded-foam sample
The graded foam sample has been analyzed with the proposed THz approach and the achieved density map has been compared with the result obtained with XRM.
To compare with the THz result, an XRM density map in the x-y plane was built using MATLAB ® through the following steps.First, each XRM scan was averaged (simple mean) in the z direction, so that a vector 1 × 3520 was obtained, corresponding to an x-coordinate line, for each y position.Then, all vectors were stacked into a 13600 × 3520 to obtain the density map in the x-y plane.Finally, the spatial resolution of this map was reduced to match the THz one, by averaging it with a square tessellation of 38 × 38 pixels, corresponding to 0.25 mm 2 (i.e., the measurement offset adopted to perform the THz survey).

Calibration curve
Table 2 reports the values of the average thickness avg and the average effective refractive index n avg for the cali- bration samples C1-C7.For each sample, the standard deviation expresses the variability among the four datasets.
The average effective refractive index values were used to build the n eff − ρ calibration curve shown in Fig. 4, obtained via linear regression, including air in the data-fitting as the origin point for this curve ( ρ = 0, n eff = 1).In the plot, the data points with their standard deviations are reported as red circles, while the data points predicted by the medium theory 29 are represented by the magenta diamonds.Note that the mixing formula given by Scheller et al. 29 has been implemented considering PP as host medium and air spherical inclusions.
The relationship between ρ and n eff resulting from the regression is given by: For data points, the absolute percentage error on the estimated density corresponding to the standard deviation is obtined as: www.nature.com/scientificreports/with ρ f taken from Table 1.The values of err are reported in Table 3. Table 4 reports the average refractive indices (n avg ) estimated for samples T1-T6 with the standard devia- tions (|σ n |) , the corresponding estimated densities (\tilde{\rho }), the absolute error on the estimated density with respect to the nominal value (|ρ f − ρ|) , the corresponding percentage error (err), the refractive indices corresponding to the nominal densities (n(ρ f )) according to Eq. ( 6) and the absolute discrepancy between the estimated n avg and n(ρ f ) (|n avg − n(ρ f )|) .The average refractive indices estimated for samples T1-T6 are reported in Fig. 4 with their standard deviations as blue crosses.
Figure 5 shows the density maps for samples T1-T6 (two for each side of each sample) achieved turning the effective refractive index n eff in each measurement point into the corresponding density as given by the Eq. ( 6).

Discussion
The relationship between the effective refractive index and the density has been built using samples C1-C7, which have been carefully manufactured in order to have a uniform density profile and morphology.The THz ToF analysis confirms that these samples are uniform, as their thickness and average effective refractive index exhibit a low standard deviation, with maximum standard deviations of 0.1 mm for the thickness and of 0.017 for the average effective refractive index, respectively.
The calibration curve is built using the estimated average effective refractive indices, revealing a linear relationship between the density and the effective refractive index in the range of interest, with a percentage error lower than 10% .It is worth noting that the absolute error on the density has not been considered when fitting the calibration curve, as this is expected to be smaller than the accuracy which can be achieved using the THz ToF analysis.
The calibration curve has been validated with samples T1-T6.The estimated n avg are in very good agreement with the calibration curve, see Fig. 4. In particular, see Table 4, for samples T2, T3, T4 and T6, the estimated average densities agree with the nominal ones, with an absolute error lower than 10 kg/m 3 and a percentage error on the estimated density of 5 %, which is consistent with the one observed for the calibration samples (C1-C7).To understand why samples T1 and T5 present larger errors, the THz retrieved density maps in Fig. 5, depicting the estimated density distribution resulting for each dataset measured for each sample, have to be observed.As expected, the density pattern of these samples is not uniform, but for T2 and T3.However, for all samples but T1 and T5, the images taken for the two sides are in very good agreement (but for a possible rotation of the sample).For T1 and T5 the two sides appear to be different, thus motivating the larger standard deviation in the estimated average refractive index and hence the larger error on the estimated density.It is worth noting that the agreement between the estimated average refractive index n avg and n(ρ f ) , i.e., the refractive index resulting from the calibration curve using the nominal average densities in Table 1, is consistent with the standard deviation for all samples.
Besides corroborating the validity of the calibration curve, Fig. 5 allows to appreciate the capability of the proposed THz ToF imaging technique of providing a spatial map of the density within PP foams.Accordingly, the technique is able to retrieve manufacturing defects of samples, which seems uniform but are found to exhibit a non-uniform distribution of the density.
The analysis of the graded-foam sample confirms the capability of the proposed THz approach of mapping a very complex profile.As can be noticed from Fig. 6, the density maps obtained via XRM and THz, respectively, are in very good qualitative agreement.The quantitative agreement expressed in terms of δ ρ also appears satis- factory, as more than 70% of the pixels have an error below 50 kg/m 3 .This is a very good agreement considering that the two images are not co-registered and the data have been acquired in different times and conditions.Moreover, it has to be taken into account: (1) the very complex distribution of the material in the direction of THz wave penetration; (2) the assumptions made for the volume averaging of both density-mapping procedures.Hence, we can conclude that the achieved result confirms that the proposed THz procedure is effective for a fast characterization of optimized foams, that is, foams with well defined 3D density maps.

Conclusion
A simple THz ToF imaging approach for estimating the complex density profile of graded PP foams has been proposed.A relationship has been experimentally obtained to link the refractive index value to the density and it has been exploited to achieve a detailed map of the density pattern of uniform and not-uniform samples from THz measurements of the effective refractive index.The proposed approach provides a 2D map, which accounts for volumetric density changes and not only for the surface ones, without requiring knowledge of the point-by point thickness of the sample under test.The effectiveness has been proved by studying the complex density distribution of an optimized 3D sample and by comparing THz results with those obtained using X-ray microscopy.The results exhibit a good qualitative and quantitative agreement.
Notably, while the case of PP foams has been considered, the approach can be immediately extended to other types of plastic foams and, being simple and effective, it is easily exploitable in laboratories involved in the development of new foaming processes, as a non-invasive tool for manufacturing process control.
Future work regards the development of a data processing procedure based on an inverse scattering approach to retrieve the actual behavior of the refractive index along the propagation path.

Figure 1 .
Figure 1.Measurement setup: (a) sketch of measurement procedure, (b) example of the measured signals referred to a metallic plane (blue curve) and a penetrable sample put on the metallic plane (black curve).

Figure 2 .
Figure 2. (a) Foam production process schematic.(b) Batch foaming apparatus: 1 to vacuum pump; 2 to gas cylinder; 3 pressure transducer output, to the data acquisition system; 4 Pt100 for measurement and control of temperature, to the data acquisition system and PID controller; 5 actuated ball valve; 6 pressure vessel; 7 electric heater input, PID controlled.(c) Mould with cylindrical cavities, 8.3 mm in height and 23.2 mm in diameter each; (d) Foam samples obtained from the process.

Figure 3 .
Figure 3.The adopted measurement set-up, with the optical sampling engine module (on the left side), whose output are connected to the THz reflection head positioned in front of the samples (right side).The sample are positioned on a flat metallic plane.

Figure 5 .
Figure 5. Density maps of T1-T6 samples.The samples with similar averaged density are shown in the same colorscale.The samples are arranged row-wise from top T1, T2, T3, T4, T5 and T6.Columns 1-2 correspond to Side 1, columns 3-4 correspond to Side 2. Axis dimensions are in mm.

Figure 6 .
Figure 6.The graded-density foam sample.(a) from top to bottom: the sample; density map achieved via THz ToF imaging and using the calibration curve; density map achieved via X-ray microscopy; absolute discrepancy between the two maps.(b) Classification of the pixels of the absolute discrepancy map.

Figure 6
Figure 6 reports the results for the graded-foam sample.The figure shows the picture of the sample, the density map achieved via THz ToF imaging and using the calibration curve ρ THz ; the density map achieved via X-ray microscopy ρ XRM and the absolute discrepancy between the two maps, δ ρ = |ρ THz − ρ XRM | .The figure also shows the histogram of the classification of the pixels of δ ρ .

Table 1 .
Foam samples parameters: sample name, density, density absolute error, volumetric fraction of the polymer.

Table 2 .
Average thickness and effective refractive index estimated for the calibration samples C1-C7.The standard deviations are also reported.

Table 3 .
Absolute percentage error for the estimated density for samples C1-C7.

Table 4 .
Test samples T1-T6: average refractive index, standard deviation, average density resulting from the calibration curve, absolute error on average estimate density, absolute percentage error corresponding to the standard deviation, average refractive index resulting from the calibration curve using the nominal density in Table1, absolute error on refractive index.